C1R, CCL2, and TNFRSF1A Genes in Coronavirus Disease-COVID-19 Pathway Serve as Novel Molecular Biomarkers of GBM Prognosis and Immune Infiltration

Glioblastoma multiforme (GBM) is a prevalent intracranial brain tumor associated with a high rate of recurrence and treatment difficulty. The prediction of novel molecular biomarkers through bioinformatics analysis may provide new clues into early detection and eventual treatment of GBM. Here, we used data from the GTEx and TCGA databases to identify 1923 differentially expressed genes (DEGs). GO and KEGG analyses indicated that DEGs were significantly enriched in immune response and coronavirus disease-COVID-19 pathways. Survival analyses revealed a significant correlation between high expression of C1R, CCL2, and TNFRSF1A in the coronavirus disease-COVID-19 pathway and the poor survival in GBM patients. Cell experiments indicated that the mRNA expression levels of C1R, CCL2, and TNFRSF1A in GBM cells were very high. Immune infiltration analysis revealed a significant difference in the proportion of immune cells in tumor and normal tissue, and the expression levels of C1R, CCL2, and TNFRSF1A were associated with immune cell infiltration of GBM. Additionally, the protein-protein interaction networks of C1R, CCL2, and TNFRSF1A involved a total of 65 nodes and 615 edges. These results suggest that C1R, CCL2, and TNFRSF1A may be used as molecular biomarkers of prognosis and immune infiltration in GBM patients in the future.


Introduction
Glioblastoma multiforme is a prevalent intracranial brain tumor, accounting for about 40%-50% of all intracranial tumors. Not only is GBM common, but it is a cancer with all sorts of frightening features and intimidating labels. First of all, GBM patients have a higher mortality rate. There are 40,000 to 80,000 new cases of GBM in China every year, with up to 30,000 deaths per year. In particular, there has been an increase in cancer patients under the age of 34, whose five-year mortality rate is the third-highest among all cancers, after pancreatic and lung cancer. Secondly, the recurrence rate of GBM is very high. Since glioblastoma is an invasive growth and invades mainly around neurons or along white matter fibers with ill-defined borders, surgery is generally difficult to remove completely. However, unresected glioblastoma is particularly prone to relapse, and relapse is often accompanied by biological malignant progression, from low-grade GBM to high-grade GBM. Third, GBM surgery is complicated. Because glioblastoma is located in the most important and complex brain of the human body, it may lead to disability or death of the patient if it is accidentally damaged in the important tissues in the brain during surgery. The difficulty and risk of surgery are very high. Glioblastoma grows rapidly, with 70% to 80% of patients lasting 3-6 months and only 10% lasting more than 1 year [1,2]. Generally, patients with glioblastoma after timely surgery can hardly live for more than two years.
With the rapid development of supercomputer technology and the continuous improvement of modern highthroughput sequencing technology, more and more largescale gene transcriptomics and related clinical databases are freely available, and it has become an increasingly vital and effective method to explore the pathological mechanisms of diseases based on bioinformatics theory. The prediction of novel molecular biomarkers through bioinformatics analysis can provide new insights into early diagnosis, survival prediction, and eventual treatment. At present, new molecular biomarkers of liver hepatocellular carcinoma, breast cancer, lung carcinoma, and adrenocortical carcinoma have been successfully predicted using largescale clinical databases and bioinformatics tools [3][4][5][6][7]. Compared with other pathological types of glioblastoma multiforme, the majority of GBM patients died from tumor recurrence, with a 5-year survival rate of less than 3%, resulting in insufficient data of patient samples in the TCGA database. Therefore, bioinformatics analysis of GBM patients has rarely been reported.
In this study, we aim to integrate information from the GTEx database and the TCGA database to identify DEGs through full-scale bioinformatics analysis. This may help to find the underlying molecular mechanisms of GBM development and may serve as biomarkers and molecular targets for the diagnosis and prognosis of GBM patients in the immediate future.

Identification and Screening of 1923
DEGs. This study was conducted according to the flow chart ( Figure 1). The GBM-GTEx and GBM-TCGA databases were used to identify a total of 1923 DEGs. Then, we screened 1118 upregulated DEGs and 805 downregulated DEGs in GBM samples compared with nontumor samples (Figure 2 and Supplementary Table S1).

GO Analysis of DEGs.
To accurately study the function of 1923 DEGs in glioblastoma, we performed GO analysis using R software. The results showed that the most enriched BP terms were modulation of chemical synaptic transmission, neurotransmitter secretion, neutrophil activation involved in immune response, synapse organization, and so on. CCs analyses revealed that the DEGs were primarily enriched in transport vesicle and collagencontaining. MF analyses revealed that the DEGs were  2 Disease Markers primarily enriched in phospholipid binding, gated channel activity, and so on (Figure 3(a)). Gene cluster analysis on BP modules with the top P value indicated that DEGs were significantly related to signal transduction processes and neutrophil activation involved in the immune response pathway (Figure 3(b)).

KEGG Pathway Analysis of DEGs.
To elucidate the role of 1923 DEGs in glioblastoma, we used the KEGG pathway analysis. The results showed that the coronavirus disease-COVID-19, prion disease, Epstein-Barr virus infection, and phagosome pathways were significantly affected with P < 0:001 (Figure 4(a)). Cluster analysis showed that 44 DEGs were significantly enriched in the coronavirus disease-COVID-19 pathway (Figure 4(b)).

Survival Analysis of DEGs.
To study the prognostic values of DEGs enriched in the pathways of coronavirus disease-COVID-19, we analyzed the overall survival of 169 GBM patients. We found that the high expression of C1R, CCL2, and TNFRSF1A genes was particularly associated with the low survival in GBM patients ( Figure 5), while other 41 DEGs enriched in the pathway of coronavirus disease-COVID-19 were not significantly connected with the overall survival of GBM patients (Supplementary Table S1). The results suggested that C1R, CCL2, and TNFRSF1A genes might be potential prognostic factors.

Expression Analysis of DEGs in GBM Patients and
Glioblastoma Cells. To detect the expression levels of C1R, CCL2, and TNFRSF1A in glioblastoma, we analyzed the      Table S2). C1R, CCL2, and TNFRSF1A mRNA levels were significantly higher in GBM samples than in normal samples (Figure 6(a) and Supplementary Table S3), which was verified by the GEPIA online tool (Supplementary Figure S1). Furthermore, the mRNA expression levels of these three genes were positively correlated (Supplementary Figure S2). We also performed experimental validation in the cell lines. Consistently, the mRNA levels of C1R, CCL2, and TNFRSF1A were strongly upregulated in U-87MG, U-251MG, and U-1118MG cells compared with those in HMC3 cells ( Figure 6(b)). The results suggested that the expression of C1R, CCL2, and TNFRSF1A might promote the development of glioblastoma.
2.6. Immune Infiltration Analysis in GBM Patients. GO analysis suggested that DEGs were significantly associated with the immune response function. To determine which specific immune cells were responsible, we performed tumor cell immune infiltration analysis. The results indicated a significant difference in the proportion of immune cells in tumor Antigen processing and presentation Graft-versus-host disease   (Figure 7(c)). In brief, these results suggested that the immune response to GBM was a complex network and was carried out in a rigorous manner.

Epstein-Barr virus infection Staphylococcus aureus infection
2.7. C1R/CCL2/TNFRSF1A Genes Were Correlated with Immune Infiltration in GBM. To determine whether C1R, CCL2, and TNFRSF1A could be used as immunotherapy targets for GBM, we studied the correlation between the mRNA levels of C1R/CCL2/TNFRSF1A and GBM immune infiltration. The C1R expression was found to be negatively correlated with B cells, CD8 + T cells, and neutrophil infiltration and positively associated with CD4 + T cells, macrophages, and dendritic cell infiltration (Figure 8(a)). The expression of CCL2 was found to be negatively correlated with CD8+ T cells and macrophage infiltration and positively associated with B cells, CD4 + T cells, neutrophil, and dendritic cell infiltration (Figure 8(b)). The TNFRSF1A expression was negatively correlated with CD8 + T cell infiltration and positively associated with B cells, CD4 + T cells, macrophages, neutrophil, and dendritic cell infiltration (Figure 8(c)). These data suggested that there might be a correlation between coronavirus disease-COVID-19 and immune infiltration of GBM. 7 Disease Markers and VCAM1, which were mainly associated with the immune system ( Figure 9 and Supplementary Table S4).

Discussion
Glioblastoma multiforme was one of the most deadly brain tumors. Surgical resection was the standard treatment for GBM, followed by chemotherapy with temozolomide (TMZ) and radiotherapy [8]. With the development of gene therapy, immunotherapy, and vaccine therapy, GBM treatment has been improved, but treatment options for controlling GBM recurrence were limited. GBM caused systemic immunosuppression and local immune dysfunction, leading to a more complex relationship between GBM and the peripheral tumor microenvironment (TME) [2,9]. TME has been found to play a role in tumor genesis, development, and migration, as well as the development of malignancy and therapeutic resistance in a growing number of studies. TME cell composition and immune cell accessibility varied greatly between GBM subtypes and patients. These factors led to the immunosuppression of GBM, which further led to the failure of immunotherapy. The identification of immune genes and immune cell types that were actively involved in TME could help clarify the general mechanism of GBM immune suppression.
Moreover, in studies that have been reported, the onset of COVID-19 caused an activation of the immune system. SARS-COV-2 bound to the host ACE2 receptor through the viral spike protein receptor binding domain (RBD) and then invaded host cells to generate an immune cascade mechanism. Mucosal-associated invariant T (MAIT) cells and γδT cells released inflammatory cytokines to respond to invasion [10]. Early response immune effector cells, including CTL and NK cells, were activated to fight the virus [11]. A number of studies have shown that the sign of COVID-19 was the cytokine storm with elevated levels of interleukin-6 (IL-6), IL-1β, complement C1r (C1R), tumor necrosis factor-alpha (TNF-α), chemokine (C-C-motif) ligand 2 (CCL2), and granulocyte-macrophage colonystimulating factor (GM-CSF) [12].
In this study, enrichment pathway analysis indicated that DEGs in glioblastoma were strongly enriched in pathways of immune response and coronavirus disease-COVID-19 (Figures 3 and 4). The high expression levels of C1R, CCL2, and TNFRSF1A in coronavirus disease-COVID-19 pathway were significantly related to the low survival in GBM patients ( Figure 5). Moreover, the mRNA expression levels of C1R, CCL2, and TNFRSF1A in glioblastoma cells or in GBM patients were found to be strongly upregulated ( Figure 6 and Supplementary Figure S1), and the mRNA expression levels of these three genes were positively associated with each other (Supplementary Figure S2). The results indicated that the expression of C1R, CCL2, and TNFRSF1A might promote the development of glioblastoma and might be used as molecular biomarkers of prognosis and immune      Disease Markers infiltration in GBM patients in the future. A recent study found the increase of CCL2 promoted carcinogenesis through esophageal mucosal inflammation [13]. Moreover, CCL2 could promote the survival and proliferation of THP-1, prostate cancer cell line PC3, renal cell carcinoma cell line 786-O and KAKI-1, and lung cancer cell line A549 [14,15]. The growth of cutaneous squamous cell carcinoma and non-small-cell lung cancer was promoted by tumor-cell-derived complement components C1r and C1s [16]. However, there are few reports about the relationship between these three genes and glioblastoma and how C1R, CCL2, and TNFRSF1A promote the occurrence and development of GBM which requires further research in the future.

NK cells activated NK cells resting T cells gamma delta T cells regulatory (Tregs) T cells follicular helper T cells CD4 memory activated T cells CD4 memory resting T cells CD4 naive T cells CD8 Plasma cells B cells memory B cells naive
Our results also found that C1R, CCL2, and TNFRSF1A gene expression was significantly associated with GBM immune infiltration. The expression of C1R, CCL2, and TNFRSF1A was closely related to B cells, CD8+ T cells, CD4+ T cells, neutrophils, macrophages, and dendritic cell infiltration (Figure 8). Recent studies found that macrophages, monocytes, NK cells, and memory T lymphocytes were activated by the CCL2-CCR2 axis, thereby stimulating the release of proinflammatory cytokines such as tumor necrosis factor-(TNF-) α, interleukin-1, and interleukin-6 [17]. Moreover, macrophages were also activated by CCL2 to secrete tissue repair factors, such as transforming growth factor (TGF)-β, platelet-derived growth factor (PDGF), and vascular endo-thelial growth factor (VEGF) [18]. TNFRSF1A encoded a member of the TNF receptor superfamily of proteins. The encoded receptor was found in soluble and membrane-bound forms and interacted with soluble and membrane-bound forms, respectively, of its ligand, TNFα [19,20]. The binding of membrane-bound TNF-α to membrane-bound receptor induced receptor trimerization and activation, playing a function in cell survival, apoptosis, and inflammation. However, it needs further study how C1R, CCL2, and TNFRSF1A affect the immune infiltration of GBM cells.
The novel coronavirus disease outbreak of 2019 (COVID-19) has emerged as the world's most serious public health threat, infecting over 1.7 million people and killing over 100,000. Recent studies have shown that cancer patients are more likely to contract COVID-19 and have higher mortality rates [12]. The novel coronavirus SARS-CoV-2 binds to human angiotensin-converting enzyme II (ACE2) through its expressed S-protein and enters the cell, while ACE2 activates the expression of CCL2, and studies have also shown that the pathogenesis of COVID-19 is closely related to the excessive release of CCL2 [21][22][23]. Our results found that the expression levels of C1R, CCL2, and TNFRSF1A in GBM cells and patients were very high (Figure 6), and 35 gene biomarkers of immune cells were significantly correlated with C1R, CCL2, and TNFRSF1A (Supplementary Table S5). We believe that the C1R, CCL2, and TNFRSF1A genes in the       Disease Markers serum (Life Technologies, USA) and 1% streptomycin and penicillin (100 μg/ml streptomycin and 100 U/ml penicillin, Life Technologies) at 37°C in 5% carbon dioxide chamber. Total RNA was extracted from cells for qRT-PCR analyses. The primers were as follows: C1R-fw: gctgcccacaccctgtatc. C1R-rv: gcccaggaacacatccaaagag. CCL2-fw: ctcgctcagccaggtaagg. CCL2-rv: actgtgggtaccacgtctgc. TNFRSF1A-fw: cctggacagaccgagtcc. TNFRSF1A-rv: ctttgtccctggtctcacc.
4.5. Immune Infiltration Analysis. The TIMER web server provides a rich resource for comprehensive analysis of immune infiltrations in different cancer types (https:// cistrome.shinyapps.io/timer/). TIMER software was used to analyze the correlation between the abundance of 6 immune infiltrates and gene expression. Moreover, we used a two-sided Wilcoxon rank sum test to detect the effect of different somatic copy number changes on tumor invasion levels.

Protein-Protein Interaction Network Construction.
STRING was a bioinformatics database that constructed the protein-protein interaction network (PPI) of DEGs based on the predicted and known PPIs. Subsequently, we analyzed the functional interactions among proteins and visualized the PPI network with Cytoscape software (version 3.7.2). 4.7. Statistical Analysis. All the data were processed and analyzed by the R software (version 4.0.2). Mann-Whitney test and t-test were used to compare the two groups of data and the optimal cut-off value generated by the R package "SURv Cutpoint" function for survival analyses. Moreover, we Figure 9: Construction of PPI network. Three key genes were marked with yellow circles. The interacted genes were marked with green circles.